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ABSTRACT 




The 3-D image reconstruction from cone-beam projections in computerized 
tomography leads naturally, in the case of radial symmetry, to the study of Abel- 
type integral equations. If the experimental information is obtained from measured 
data, on a discrete set of points, special methods are needed in order to restore 
continuity with respect to the data. A new combined Regularized-Adjoint-Conjugate 
Gradient algorithm (introduced in this work), together with two different 
implementations of the Mollification Method (one based on a data filtering 
technique and the other on the mollification of the kernel function) and a 
regularization by truncation method (initially proposed for 2-D ray sample schemes 
and more recently extended to 3-D cone-beam image reconstruction) are extensively 
tested and compared for accuracy and numerical stability as functions of the level 
of noise in the data. 


1. INTRODUCTION. 

The difficult problem of determining the structure of an object from its 3-D 
cone-beam data projections is currently receiving considerable attention (see B. D. 
Smith, Ref [16]). When the object is known to be radially symmetric, its structure 
can be determined by using the inverse Abel transform. If the object does not have 
radial symmetry, it can be reconstructed, in principle, by using the inverse Radon 
transf orm. 

Abel's integral equation can be written as 


f(x) « S g(s) (x - s)" 1/2 ds, 0 i x s 1, 
o 


( 1 ) 


where the function f(x) is the data function and g(s) is the unknown function. The 
exact solution is given by 


1 X 

g(x) = - S f’(s) (x - s)’ 1/2 ds, 0 £ x £ 1, 
n o 


( 2 ) 


provided the derivative exists and f (0) = 0. (See R. Gorenflo and S. Vessella, Ref 

[ 6 ]). 

It is well-known (References [1], 12], [4] and [6]) that Abel’s integral equation 
is somewhat ill-posed, that is, small errors in the data f(x) might cause large 
errors in the computed solution g(x). Consequently, the direct use of formula (2) 
is very limited and special methods are needed. 
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r This paper has two main purposes. First, we present and briefly analyze a new 

stable method for the numerical solution of Abel’s integral equation, Method I, by 
weakly coupling the original problem with its adjoint formulation obtaining a 
regularized system of linear equations which is then successfully solved by the 
conjugate gradient method. Second, we test and compare the numerical stability and 
the accuracy of Method I and three other known algorithms on several benchmark 
examples as a function of the amount of noise in the data. 

Method II in this paper (see D. A. Murio, Ref [12]), is obtained by initially 
filtering the noisy data by discrete convolution with a suitable averaging kernel 

instead of mollifying the kernel function in equation (2), Method III, as required 
by K. Miller (Refs [10] and [11]) in his reconstruction algorithm for 2-D ray- 
sampling schemes. Method IV has been implemented by D. A. Murio, D. Hinestroza and 
C. E. Mejia (Ref [13]) based on a regularization by truncation technique initially 
proposed by B. K. P. Horn (Ref [9]) and recently extended to 3-D image 

reconstruction methods from cone-beam projections by B. D. Smith (Ref [15]). 

In Section 2 we introduce the new Method I, analyze the consistency and stability 
properties of the algorithm and obtain an upper bound for the error. In Section 3, 
we describe the other procedures and discuss in detail the numerical implementation 
of all the methods involved. Section 4 is devoted to the numerical testing of the 
four algorithms and the presentation of several useful comparisons involving 
Methods I, II, III and IV. Some conclusions are included in Section 5. 


2. REGULARIZED-ADJOINT-CONJUGATE GRADIENT METHOD. (Method I). 

In a more abstract setting, equation (1) can be written as 

Ag = f, 

where A represents the Abel integral operator. For suitable f unctions h and q, the 
adjoint operator A* is defined by 


A h(x) = q(x) a J h(s) (s-x)' 1/z ds, 


0 s x s 1, 


and it is clear that the homogeneous equation A*h = 0 has the unique solution h(x) 
s 0, 0 s x s 1. Hence, as a direct consequence of Fredholm alternative (see P. R. 
Garabedian, Ref [5]), solving the singular equation Ag = f for smooth but 
otherwise arbitrary data functions f satisfying f(0) = 0, is equivalent to solve 
the uncoupled system of linear integral equations 

_ Ag = f 

A* h = 0. (3) 

In order to help stabilize the inverse problem, we propose to solve, instead of 
(3), the weakly coupled system of equations 


r 

Au - av = f 

" A*v + au = 0, 0 < a « 1, (4) 

V 

by successive approximations. This system is equivalent to 


L 
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r 

r 

Au - av = f 

a/3A*v + a 2 /3u - u + u = 0, 


where j3 is any nonzero real number to be determined. We elect fi to depend on the 
iteration and rewrite the previous system as 


Au n - av n = f 


a M* v n + a 2 Pn u n “ u n + u n+l = °- 


to obtain 


_ av n = Au n " f 

u n +i = u n - /3 n (a 2 u n + A*(av n )l, n = 0,1,2,..., (5) 

w 

u 0 arbitrary, usually 0. 

Remarks : 

1. Each iteration in (5) involves the solution of two ’’direct" problems: one 
corresponding to the original operator, Au n , and the other associated with the 

adjoint operator, A*(av n ). 

2. Elimination of v in system (4) leads to the set of normal equations, with I 
indicating the identity operator, 


(A*A + a z I)u = A*f , (6) 

which characterizes the minimum of the zero order Tikhonov functional (see C. W. 
Groetsch, Ref [7]) 

J(u) = - ( 1 1 Au - f 1 1 2 + a 2 1 1 u 1 1 2 ). (7) 

2 

3. The gradient of the functional (7) is given by 

VJ(u) = a 2 u + A*(Au - f) 

and it is easily computed if the solution of the adjoint problem is known. In fact, 

taking into consideration (4), we can write VJ(u) = a 2 u + A*(av), and for each 
iteration we get 

VJ(u n ) = a 2 u n + A*(av n ). (8) 

These considerations allow us to choose |3 nf for each n, in such a manner that 
system (5) can now be solved by the Conjugate Gradient Method (W. M. Patterson, Ref 
[14]). 

The complete abstract algorithm, after introducing the notations 

l 

(f,g) = J f (x)g(x) dx and II T II = (f,f) 1/2 , corresponding to the inner product and 
o 
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r norm respectively of square integrable functions on the interval [0,1] is as 
follows: 


0 ) 

1) 

2 ) 

3) 

4) 

5) 

6 ) 

1 ’) 

2 ’) 

3’) 

4’) 


4”) 


5’) 

6 ') 


For n = 0, 

Set u 0 = 0 and choose a > 0. 

Compute Au 0 , i.e., solve the original direct problem. 
Compute the residual av 0 = Au 0 - f. 

Compute A*(av 0 ), i.e., solve the direct adjoint problem. 
Evaluate the gradient d 0 = VJ(u 0 ) using formula (8). 

1 1 d 0 1 1 2 

Set r 0 = . 

a 2 1 1 d 0 1 1 2 + 1 1 Ad 0 1 1 2 
Update: Uj = u 0 - r 0 d 0 . 

For n = 1,2 

Solve the original direct problem Au n . 

Compute the residual av n = Au n - f. 

Solve the direct adjoint problem A*(av n ). 

Evaluate the gradient VJ(u n ) using formula (8). 

II VJ(u n ) II 2 

Compute d n = VJ(un) + rf n _,. 

II VJ(u n _j) II 2 
(VJ(u n ),d n ) 

Set r n = . 

a 2 1 1 d n 1 1 2 + 1 1 Ad n 1 1 2 
Update: u n+1 = u n - r n d„. 


Stability of Method I. 

We consider now the more realistic situation when instead of the exact data 
function f, we only know some noisy data function f e satisfying 

1 1 f - f e 1 1 s e. 

In this section the unique solution of system (4) will be denoted by ug to 

emphasize its dependency on the regularization parameter a and the level of noise 
in the data c. Assuming that the ideal problem (1) for errorless data f has the 
unique solution g = A*‘f, since ug satisfies equation (6), with f replaced by f E , 
it follows from well-known estimates in the theory of Tikhonov regularization that 

1 1 g - u a 1 1 = C 0 a l/2 , and 1 1 u a - ug 1 1 s e cr‘ /2 , 

for some constant C 0 > 0, independent of a; u a denotes the regularized solution 
when c = 0. 

Combining these estimates, we obtain the error upper bound 


L 



243 

Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 

riCTDF.S-Iin Editor: G.S. Dulikravich. Washington D.C.. October 212L 199 L z 

r 

1 1 g - u G I I s C 0 a 1/2 + e a 1/2 

and choosing a = CjC for some constant C : > 0, it follows that 

II g - ug II s (C 0 + Cj 1 ) a 1/2 (9) 

which shows that, theoretically, as the quality of the data becomes better and 
better (c -> 0), we get convergence with rate a 1/2 . See C. W. Groetsch (Ref [7]) for 
details. 

The convergence of the sequence of iterates ug <n from system (5), with 0 n as 

discussed above, to the unique solution u£ of the canonical equations (6) as n -» oo 
is well documented, for instance, in the work of C. W. Groetsch, J. T. King and D. 
A. Murio (Ref [8]) and will not be pursued further here. 

The finite dimensional version of the combined Regularized-Adjoint-Conjugate 
Gradient algorithm will be discussed in the next Section. 

3. METHODS II, III AND IV. NUMERICAL IMPLEMENTATIONS. 

Method II in this paper is based on attempting to reconstruct a mollified version 
of the solution g in equation (2). After introducing the 6-mollifier 

p Jx) = - 7 T 1/2 exp[-x 2 /5 2 ] (10) 

6 6 

of "blurring radius" 5 and extending the data function f G to the interval 
[-35,1+36] in such a way that it decays smoothly to zero on [1,1+35] and it is zero 
on [-36,0], an approximate solution is defined by 

gf(x) = - JMp * f c r(s) (x - s)- 1/2 ds, o S X S 1. (11) 

6 71 0 6 


Here, 

(p * f G )’(x) = J°° — IpJx - s) f c (s)] ds « f 3 — [p Ax - s) f c (s)] ds, 

O —oo dx 6 X-36 dx d 

showing that the main idea of the method consists on replacing the noisy data 
function f c by the filtered data function p^* f G . It is important to notice that 

the radius of mollification, 5, can be uniquely and automatically determined as a 
function of the amount of noise in the data, c, based in the fact that there is a 
unique value of the regularizing parameter 6 for which 

II p * f e - f G II = c. (12) 

o 

Under very mild conditions, i.e., if f e is continuous and if the second 
derivative of the errorless data function f is uniformly bounded by M 2 in the 


244 


Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 

H CIDE S-rii ). Edit o r; O. S . P uii k ravj gh . W a shin gton P .C ., QsMzi21 ^ l ML 

sample interval (0,1), the following error estimate holds 
e 6 

I I g g - g 1 1 5 - (6 M 2 + c/5). ( 13 ) 

The complete abstract algorithm is as follows: 

1) Automatically determine the unique radius of mollification 6 as a function of 
the level of noise c. 

2) Smoothly extend the noisy data function f e to 1-35,1+35], 

3) Compute the derivative of the filtered data function p * f e . 

6 

4) Compute g G using equation (11). 

For more details and further discussions, the reader should consult D A Murio 
Ref [12]. 

Method III is based on the Mollification Method as originally proposed by K. 
Miller (Refs (10] and [11]) for 2-D ray-sampling reconstruction geometries. First 
we notice that the exact formula (2) can be written 

1 

g(x) = - (k * f’)(x), 0 a x s l, 

71 

where k(t) = r 1/z represents the kernel function. The mollification of the last 
equation with the averaging kernel defined in (10) gives 

1 

(p * g)(x) = - (p * k * f’)(x). 

° n ° 

In Method II, we associated the right-hand side of this equation as k* (p g * f’) = 
k* (p g * f )’ ; for Miller’s idea we associate as (p g * k) * f’= (p g * k)’* f and obtain 
the approximate reconstruction solution 

1 X 

= - S (p~* k)’(x - s) f G (s) ds, 0 a x a 1. (14) 

° TT 0 ° 

Mathematically, formulae (11), for Method II, and (14), for Method III, are 
identical. Consequently, the theoretical error bound (13) derived for Method II 
also applies for Method III. 

The complete abstract algorithm for Method III is given by: 

1) Choose 5 > 0. 

2) Compute the mollified kernel p * k. 

O 

3) Evaluate the derivative of the mollified kernel p * k. 

o 

4) Compute g| using equation (14). 

Remarks: 

1. In Method III, the mollified kernel is computed only once and is used repeatedly 
for different data functions. 


L 


J 
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| 2. Method II requires a filtering of each data function and the corresponding 

parameter is automatically selected according to the quality of the measured data. 

3. The selection of the mollification parameter in Method III requires further 
consideration. 

Method IV is based on a reconstruction technique initially proposed by B. K. P. 
Horn (Ref [9]) for arbitrary 2-D ray schemes and more recently extended to 3-D 
image reconstruction methods from cone-beam projections by B. D. Smith (Ref [15]). 
Integrating by parts equation (2), we obtain the equivalent expression 

g(x ) = - lim / TT 3/2 S f(s) ds S *f(s)(x-s)~ 3/2 ds V 0 - x ^ 1. 

Tt V- x-y 2 o J 

The approximate inverse Abel transform is now obtained by eliminating the limit 
procedure in the last expression, i.e., 

gC (x) = L / y-3/2 J x fC( s ) ds - - rV(s) (x - s)- 3/z ds 0 ± X s 1. (15) 

y TT ^ x-y 2 0 ' 

By requiring the second derivative of the errorless data function f and the 
measured data function f c to be continuous, we obtain the following error estimate 

5 y l/2 2 

1 1 g e - g 1 1 s - Mj + - e r 1/2 + 0 (* 3/2 ), 06) 

y 2 7T TT 

where Mj is a uniform bound for f* on the interval (0,1). For a proof of this 
assertion and a complete analysis of Method IV, see D. A. Murio, D. Hinestroza and 
C. E. Mejia (Ref [13]). 

The complete abstract algorithm for Method IV is reduced to 

1) Choose * > 0. 

2) Compute g c using formula (15). 

7f 

Remark : 

The error estimates (9), (13) and (16) show that all the methods are consistent 
and stable with respect to perturbations in the data, in the L 2 norm, for a fixed 

choice of the several regularization parameters a, 6 or y. 

Numerical Implementations . 

Since in practice only a discrete set of data points is generally available, we 

assume that the data function f c is a discrete function measured at equally spaced 
sample points on the interval [0,1]. For h > 0 and Nh = 1, we let Xj = jh and 

denote f c (Xj) = f j, j = 0,1, with f£ = 0. 

Method I: 

Discretization leads to a finite dimensional version of the combined 
Regularized-Adjoint-Conjugate Gradient algorithm of Section 2. The operators A and 

A* are represented now by a matrix A and its transpose A T , respectively. The 

L 
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r . * ] 

approximate discrete solution u mfCt , obtained after m iterations, the gradient 

VJ(u m ), d m , r m , u 0 and the residual av m are now N-dimensional real vectors. From 

equation (1), a simple discretization gives the lower triangular system of linear 
equations 

h a j+l-l^ U m,a^i = fj * 

1=1 J 

where 

aj = ( jh)~ 1/z , j = 1,2,. ...N. 

indicates the ( j-1) subdiagonal of the N x N matrix A. 

The discrete algorithm for the Conjugate Gradient method (see P. G. Ciarlet, Ref 
[3]) follows exactly the steps described previously in Section 2, and we only have 
to add the necessary stopping criteria, given by 

I I u m,a " u m-i,a ^ 2 - TOL 1 1 Uj^^. 1 1 2 , 

where TOL is a small positive tolerance parameter entered by the user and 

II f It 2 = { - I N If jl 2 

v N j=o J 

is the discrete 1 2 norm on [0,1], 

Method II: 

To numerically approximate g^(x), a quadrature formula for the convolution 

equation (11) is required. The objective is to introduce a simple approximation and 
avoid any artificial smoothing in the process. 

Given Xj, j = 0,1 N, we define 

q E (x) = £ J f^,(x), 0 < x < x., 

1=0 3 


y/2 


(17) 


a piecewise constant interpolation of f c (x) at the grid points Xj. Here, 


4>o (x) = { 


1, Oi x £h/2 
0, otherwise ’ 



Xj-h/2 ^ x 
otherwise 


x 


j 


and 


4>i(x) = 


1, X(-h/2 £ x £ x,+h/2 
0, otherwise 


i = 1,2,..., j-1. 


The computational algorithm is as follows: 

After smoothly extending the discrete data function to any interval of interest 
containing the sample interval [0,1], we determine the radius of mollification <5 as 
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[~ a function of the amount of noise in the data e by solving the discrete version of 
equation (12) using the bisection method. Next, we substitute f c by its 
interpolation q e and compute the approximation to p^ * f c given by the discrete 

convolution 

(p g « q E )(x) = £ f E (p g « 0 k )(xj) = £ f E m^. 

where the weights m^ are evaluated exactly. A discrete version of the derivative of 

the discrete filtered data function is obtained using centered finite differences. 

Finally, the discrete approximation to g^ is calculated by discretely convolving 

the computed derivative approximation against the sampled data function (see 
equation (ID). For a detailed analysis of this algorithm, the reader is referred 
to D. A. Murio (Ref [12]). 

Method III : 

The convolution p- * k requires an extension of the singular kernel k for values 
o 

of x less or equal to zero. In our implementation we use the following symmetric 
extension: 


k(0) = 2h“ 1/2 , k(-x) = k(x), x > 0. 

The discrete approximation is now straightforward: 

With Sj = (p^ * k)(Xj), j = 0,1, ...,N, the discrete convolution formula 

corresponding to equation (14) is 

g£ >5 (0) = 0. 

gU(x ,) « I E J Sj_k( f k + 1 - fk-,)/2. j = 1.2....N-1 

J TT k=l 

and 

g£,s(D = g + - s 0 (fS - fS-iJ. 

71 

where is the approximate inverse Abel transform at the grid points. 

Method IV : 

In this case, we first construct a piecewise linear interpolation of f G (x) at 
the grid points x jf given by 

q e (x) = J] J f j 0i(x), 0 < x ± Xj, 

j 


where the functions ^(x), i = 0 , 1 ,..., N are given by 

L 
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r 


a fvi — I 1-x/h, 0 5 x 5 h 
0 l 0 otherwise, 


0,(x) = / 1+(x " x J )/h - Xj.^x^xj 
3 l 0 otherwise, 


and 


f 1+(x-x i) /h - X 1-1 s X 5 X 
•Pi (x) ~ \ l-(x-Xj)/h, x ( 5 x 5 Xj + 
0 otherwise, 


i = 1.2 j. 


We notice that the approximate solution g|(x) of formula (15) can also be 
written as 


gf(x) = - (H * fC)(x), 0 5X51, 

TT 


where the kernel is defined by 

Hy(t) 


r' 3/2 . 0 5 t < y 

4 1 

t' 3/2 , ? < t. 

2 


(18) 


The quadrature formula for equation (18) is obtained by directly convolving the 
kernel function with q e as indicated below. Thus, the computed solution at the 

grid points is given by 

gf.h( x j> - - q e )(x,) = - l J f® bf(x,), 

n TT i=0 J 

where the weights 

bfUj) = S J H^(Xj - s) ^j(s) ds 

are evaluated exactly for i = 0,1,. ..,j. The readers interested in further details 
should consult D. A. Murio, D. Hinestroza and C. E. Mejia (Ref [13]). 

4. NUMERICAL RESULTS AND COMPARISON. 

In this section we describe the tests that have been implemented in order to 

compare the performance of the methods introduced in previous sections. 

We tested the methods on three examples. In all of them, the exact data 
function is denoted f(x) and the noisy data function f c (x) is obtained by adding an 

c random error to f(x), that is, f c ( Xj ) = f( Xj ) + c<r } , where Xj = jh, j = 

0,1,. ..,N; Nh = 1 and cr } is a uniform random variable with values in [-1,1] such 
that 

max I f e (x,) - f(x.) I ^ c. 

0^j<N' j J 1 




L 


The exact inverse Abel transform is denoted g(x) and its approximation given by 
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any of the methods is denoted g£ h (x), where p represents the regularization 
parameter of the particular method. 

Example 1: 

As a first example we consider the data function f(x) = x with exact inverse 

2 

Abel transform g(x) = - x 1/2 . This data function satisfies all the necessary 

7T 

hypotheses for convergence estimates of Sections 2 and 3. 

Example 2: 

The data function 


f(x) = 


2x 2 

l-2(l-x) 2 ! 


0^x< 1/2 
1/2 sx^l, 


is only once continuously differentiable on [0,1], partially violating the required 
conditions for the theoretical error analysis of Sections 2 and 3. In this example, 
the exact inverse Abel transform is given by 


(16/3n)x 3/2 , Osx 1/2 

( 16/3n )x 3/2 + (16/3tt) ( x-l/2) 3/2 - (8/Tr)(x-i/2) l/2 (2x-l), 
w 1/2 <X— 1 . 


Example 3: 

The data function is defined as follows: 


f(x) = 


r 

o. 

- 2(x-0.2) 1/z , 
2(x-0.2) 1/z - 


2(x-0.6) 1/z 


0 i x <0.2, 
0.2 £ x s 0.6, 
0.6 < x si. 


Its first derivative is not continuous on 10,1], strongly violating the 
necessary hypotheses for the convergence estimates of Sections 2 and 3. The exact 
inverse Abel transform is given by 


g(x) 


1, 0.2 s x s 0.6, 

■< 

0, otherwise. 


The four methods were tested for three different values of N, N = 200, 500 and 
1000, three different values of e, e = 0.0, 0.005 and 0.01, and several values of 
the corresponding regularization parameters. The algorithms were extensively used 
and we numerically determined appropriate values for the regularization parameters 
for each method, except for Method II where the radius of mollification was 
selected automatically. These quasi-optimal parameter values are used in the tables 
,and figures below. 
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Different values of e provide a crucial test for stability. Tables 1, 2 and 
illustrate this point. The error norms in the tables are computed as II g - g e || 

according to definition (17). In the tables, each row corresponds to one of the 
methods with a fixed regularization parameter, and shows the change in the error 
norm due to changes in the level of noise in the data. The presented numerical 
results indicate stability. The columns in the tables allow us to compare the 
performance of the methods under similar conditions. 

Figures 1 to 4 show the reconstructions of the step function of Example 3 
provided by the four methods for the same number of sample data points, N = 500, 
the same noise level, c = 0.01, and quasi-optimal regularization parameters. The 
qualitative behavior is quite good taken into consideration the high amount of 
noise in the data. 


1 


5. CONCLUSIONS 

The following are some conclusions based on the implementations of the methods 
presented in this paper: 

Consistency and stability of the four methods is clearly confirmed throughout 
experimentation and very weak dependency on the parameter N is observed. 

Method II provides an automatic mechanism to select the radius of mollification 
as a function of the level of noise in the data. Furthermore, as a consequence of 
the stability of the four methods, it is easy to find, by numerical 

experimentation, lower and upper bounds for quasi-optimal regularization parameters. 

An advantage of method III over method II is that the mollification of the 
kernel is computed only once and can be used for different data functions. Methods 
II applies mollification to each data set. 

All the results are very competitive. However, mollification solutions are 
slightly better in terms of accuracy and method IV, the easiest to implement, seems 
to be more sensitive to perturbations in the data. 
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Method 

Parameter 

c = 0.0 

c = 0.005 

c = 0.01 

I 

a 

= 0.08 

0.0279 

0.0294 

0.0359 

II 

6 

= 0.008 

0.0000 

0.0048 

0.0096 

III 

6 

= 0.008 

0.0005 

0.0137 

0.0274 

IV 

7 

= 0.004 

0.0302 

0.0315 

0.0349 


Table L Error Norms as functions of e 
in Example 1 with N = 500 


Method 

Parameter 

c = 0.0 

c = 0.005 

c = 0.01 

I 

a = 0.08 

0.0275 

0.0293 

0.0365 

II 

6 = 0.008 

0.0001 

0.0048 

0.0096 

III 

6 = 0.008 

0.0005 

0.0136 

0.0273 

IV 

7 = 0.001 

0.0174 

0.0263 

0.0431 


Table 2. Error Norms as functions of c 
in Example 2 with N = 500 


Method 

Parameter 

e = 0.0 

c = 0.005 

e = 0.01 

mm 

a = 0.08 

0.0615 

0.0618 

0.0641 

K^l 

6 = 0.008 

0.0052 

0.0052 

0.0053 

Bpg*, 

6 = 0.008 

0.0295 

0.0330 

0.0411 


r = 0.001 

0.0648 

0.0678 

0.0760 


Table 3. Error Norms as functions of c 
in Example 3 with N = 500 
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Fig. 1 Reconstruction with Method I 
C = 0.01, a = 0.08, N = 500 


Fig. 2 Reconstruction with Method II 
€ = 0.01, 6 = 0.008, N = 500 
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Fig. 3 Reconstruction with Method HI 
c = 0.01, 6 = 0.008, N - 500 


Fig. 4 Reconstruction with Method IV 
C = 0.001, 7 = 0.001, N = 500 


